* Creates Figure 1 and Figure A3

clear all

use "$savedata/masterdata.dta", replace

keep if sample25==1
gen vol = vol25

reghdfe survive30 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe30, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe30 = docfe30*signal_2step

drop residual* sigma* signal*

reghdfe survive365 c.prevyear_cost, absorb(i.sex##i.derv_age i.black i.mixed i.chinese i.asian i.race_miss i.ynch* i.prevyear_stroke i.di1 i.di2 i.di3 i.di4 i.di5 i.shock i.arythmia i.arthero i.arrest i.dow##i.admidate_mont##i.finyear i.hyid, savefe) keepsingleton
predict residual, residuals

xtreg residual c.std_ami3 c.std_nonami3, fe i(doctor_id)
predict docfe365, u
gen sigma_u = e(sigma_u)
gen sigma_e = e(sigma_e)

gen signal_2step = ((sigma_u^2) / ((sigma_u^2) + ((sigma_e^2)/vol)))
gen adj_docfe365 = docfe365*signal_2step


collapse (mean) *docfe* vol, by(doctor_id)

*Creates Figure A3
scatter adj_docfe365 adj_docfe30, mcolor(gs3) || lfit adj_docfe365 adj_docfe30, ytitle(Estimated doctor quality (365 days)) xtitle(Estimated doctor quality (30 days)) graphregion(color(white)) ///
legend(off) 
graph export "$results/figA3.pdf", as(pdf) replace

* Creates Figure 1, panel A (3o day survival)
twoway kdensity docfe30, lcolor(black) || kdensity adj_docfe30, lcolor(black) lpattern(dash)  xtitle(Estimated doctor FE (30-days)) graphregion(color(white)) ///
legend(off) xlab(-0.2(0.1)0.2) ytitle("Density")
graph export "$results/fig1a.pdf", as(pdf) replace

* Creates Figure 1, panel B (365 day survival)
twoway kdensity docfe365, lcolor(black) || kdensity adj_docfe365, lcolor(black) lpattern(dash)  xtitle(Estimated doctor FE (365-days)) graphregion(color(white)) /// 
legend(label(1 "Unadjusted") label (2 "Adjusted")) xlab(-0.2(0.1)0.2) ytitle("Density")
graph export "$results/fig1b.pdf", as(pdf) replace

